Quality Control
jeff karnsomprot
2025-06-05
1 Objective
The goal of this session is to get a clean macula densa dataset to work with.
2 Load Packages
options(future.globals.maxSize = 74 * 1024^3) # 55 GB
getOption("future.globals.maxSize") #59055800320## [1] 794568949764 Finding Non-MD Cells
SO1 <- SCTransform(SO1) %>%
RunPCA() %>%
FindNeighbors(dims = 1:40) %>%
FindClusters(resolution = 2) %>%
RunUMAP(dims = 1:40)## Running SCTransform on assay: RNA## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 23047 by 11777## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07## Found 280 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 23047 genes## Computing corrected count matrix for 23047 genes## Calculating gene attributes## Wall clock passed: Time difference of 26.12352 secs## Determine variable features## Centering data matrix## Place corrected count matrix in counts slot## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Set default assay to SCT## PC_ 1
## Positive: Klf2, Tm4sf1, Ifitm3, Vim, Crip1, Tmsb4x, S100a6, Ly6c1, Fbln5, Slc9a3r2
## Ptprb, Mgp, Stmn2, Cldn5, Ltbp4, Srgn, Icam2, Podxl, Id1, Egfl7
## Pecam1, Epas1, Pi16, Jag1, Fxyd5, Fabp4, Clec14a, Gja4, Cavin1, Eng
## Negative: Umod, Egf, Klk1, Atp1b1, Wfdc2, Atp1a1, Sfrp1, mt-Nd4l, Ldhb, mt-Nd5
## Slc12a1, Sostdc1, Gm13340, mt-Co1, Ppp1r1a, Gm28437, Mt1, Tmem213, Tmem52b, Gm29216
## Kng2, Car15, Cd63, Cox8a, Spp1, Fxyd2, Gm28661, Sdc4, mt-Cytb, Gm10925
## PC_ 2
## Positive: Umod, Egf, Tmem52b, Sult1d1, Fabp3, Cox6c, Krt7, Foxq1, Sostdc1, Prdx5
## Klk1, Wfdc15b, Wfdc2, Slc25a5, Ly6a, Ckb, Atp5g1, Cox5b, Cldn19, Cox4i1
## Ndufa4, Cox8a, Atp1a1, Chchd10, Ggt1, Gm47708, Gm44120, Igfbp7, Gm10053, Pcolce
## Negative: Pappa2, Malat1, Zfand5, Aard, Gm37376, Robo2, Wwc2, Pde10a, Neat1, Itga4
## Sdc4, Nadk2, Nos1, Irx1, Ramp3, Sgms2, Col4a3, Mir6236, Ptgs2, Col4a4
## Tmem158, Etnk1, Bmp3, Ranbp3l, Itprid2, Sfrp1, Cdkn1c, Sec14l1, Camk2d, Aebp1
## PC_ 3
## Positive: Fos, Junb, Jun, Egr1, Btg2, Hspa1b, Hspa1a, Atf3, Zfp36, Ier2
## Fosb, Socs3, Klf6, Jund, Klf2, Dnajb1, Dusp1, Gadd45b, H3f3b, Tsc22d1
## Gadd45g, Cebpd, Ubc, Actb, Rhob, Mt1, H1f2, H2bc4, Mt2, Ppp1r15a
## Negative: Egf, Umod, Mir6236, CT010467.1, Slc12a1, Pappa2, mt-Nd5, mt-Nd4l, Ly6c1, Gm37376
## Ptprb, Ltbp4, Pecam1, Stmn2, mt-Rnr1, Fbln5, Cldn5, Eng, Icam2, Egfl7
## mt-Co1, Mgp, Epas1, Podxl, Heg1, Ramp2, Slc9a3r2, Sox18, Fbln2, Pi16
## PC_ 4
## Positive: Gm22133, Fth1, Ldhb, Ubb, Ftl1, Cd63, Car15, Prdx1, Mt1, Mpc2
## Eif1, Cd9, Fxyd2, Mgst1, Tmem213, Rps24, Bsg, Clu, Itm2b, Rpl7
## Rps5, Mdh1, Rps27a, Rpl9, Aldoa, Rpl11, Rps7, Wfdc2, Ramp3, Selenow
## Negative: Malat1, Mir6236, Gm37376, Egf, CT010467.1, Umod, Fos, Neat1, Nme7, Tmem52b
## Junb, Jun, Gm24447, mt-Nd5, mt-Nd6, mt-Rnr1, Egr1, mt-Nd4l, Kcnq1ot1, Slc12a1
## Atf3, Fosb, Zfp36, Btg2, Lars2, Klf6, Socs3, Dst, Slc5a3, Atp1b1
## PC_ 5
## Positive: Umod, Egf, Mt1, Ly6a, Hspa1a, Hspa1b, Jun, Ppp1r1a, Cldn10, Sostdc1
## Atp1a1, Kng2, Mt2, Gm22133, Cd63, Sfrp1, Ckb, Car15, Klf2, Fos
## Junb, Slc12a1, Wfdc15b, Pth1r, Klf6, Wfdc2, Cldn19, Socs3, H3f3b, Ggt1
## Negative: Rhcg, S100g, Aqp3, Spink8, Calb1, Pgam2, Hsd11b2, Scnn1b, Aqp2, Rhbg
## Fxyd4, Defb1, Cav2, Kcne1, Cryab, Slc8a1, Scnn1g, Apela, Iqgap2, Tmem52b
## Tmem45b, Tmsb4x, Tbx2, Cav1, Npnt, Gata3, Tmem229a, Car2, Gm43305, Igfbp7## Computing nearest neighbor graph## Computing SNN## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11777
## Number of edges: 396350
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7209
## Number of communities: 29
## Elapsed time: 1 seconds## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 21:38:40 UMAP embedding parameters a = 0.9922 b = 1.112## 21:38:40 Read 11777 rows and found 40 numeric columns## 21:38:40 Using Annoy for neighbor search, n_neighbors = 30## 21:38:40 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 21:38:40 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpumjlXP/file671a2a894134
## 21:38:40 Searching Annoy index using 1 thread, search_k = 3000
## 21:38:42 Annoy recall = 100%
## 21:38:42 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 21:38:43 Initializing from normalized Laplacian + noise (using RSpectra)
## 21:38:43 Commencing optimization for 200 epochs, with 519266 positive edges
## 21:38:43 Using rng type: pcg
## 21:38:45 Optimization finished
# Viewing with Known Kidney Celltype Markers
markers.to.plot1 <- c("Lrp2", # PT
"Slc5a12", # PT-S1
"Slc13a3", # PT-S2
"Slc16a9", # PT-S3
"Havcr1", # Injured PT
"Epha7", # dTL
"Cryab", # dTL
"Cdh13", # dTL1
"Slc14a2", # dTL2
"Slc12a1", # TAL
"Umod", # TAL, DCT1
"Egf", # TAL, DCT1,
"Cldn10", # TAL
"Cldn16", # TAL
"Nos1", # MD
"Slc12a3", # DCT
"Pvalb", # DCT1
"Slc8a1", # DCT2, CNT
"Aqp2", # PC
"Slc4a1", # IC-A
"Slc26a4", # IC-B
"Nphs1", # Podo
"Ncam1", # PEC
"Flt1", # Endo
"Emcn", # Glom Endo
"Kdr", # Capillary Endo
"Pdgfrb", # Perivascular
"Pdgfra", # Fib
"Piezo2", # Mesangial
"Acta2", # Mural
"Ptprc", # Immune
"Cd74", # Macrophage
"Skap1", # B/T Cells
"Upk1b", # Uro
"Top2a" # Proliferation
)
DotPlot(SO1,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b
5 Reclustering
# Reclustering after Subsetting Non-MD Cells
SO2 <- SCTransform(SO2) %>%
RunPCA() %>%
FindNeighbors(dims = 1:36) %>%
FindClusters(resolution = 2) %>%
RunUMAP(dims = 1:36)## Running SCTransform on assay: RNA## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 14607 by 11455## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07## Found 274 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 14607 genes## Computing corrected count matrix for 14607 genes## Calculating gene attributes## Wall clock passed: Time difference of 16.34481 secs## Determine variable features## Centering data matrix## Place corrected count matrix in counts slot## Warning: Different cells and/or features from existing assay SCT## Set default assay to SCT## PC_ 1
## Positive: Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Klk1, Prdx5, Cox6c, Foxq1
## Wfdc15b, Ly6a, Wfdc2, Slc25a5, Krt7, Cox5b, Ckb, Atp5g1, Cldn19, Cox4i1
## Atp1a1, Ndufa4, Ggt1, Cox8a, Chchd10, Atp5b, Cox6b1, Uqcrq, Uqcr11, Cox7a2
## Negative: Pappa2, Zfand5, Aard, Robo2, Wwc2, Pde10a, Itga4, Nadk2, Mir6236, Nos1
## Neat1, S100g, Ramp3, Col4a4, Sgms2, Col4a3, Ptgs2, Irx1, Tmem158, Itprid2
## Ranbp3l, Bmp3, Camk2d, Sdc4, Cdkn1c, Etnk1, Zbtb20, Srrm2, Acsl4, Hnrnpa2b1
## PC_ 2
## Positive: Fos, Junb, Jun, Hspa1b, Hspa1a, Btg2, Egr1, Zfp36, Atf3, Ier2
## Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Rhob, Tsc22d1, Gadd45g
## H3f3b, Gadd45b, Ubc, Cebpd, Actb, H1f2, Mt1, Mt2, H2bc4, Klf4
## Negative: Pappa2, Mir6236, CT010467.1, Egf, Slc12a1, Umod, Etnk1, Sfrp1, Wnk1, Atp1b1
## Robo2, Hsp90b1, Sec14l1, Sptbn1, Col4a4, Nme7, Pde10a, Itprid2, Zbtb20, Mal
## Wwc2, Col4a3, App, Dst, Ly6a, Nfe2l1, Trim2, Tfcp2l1, Pou3f3, Bmp3
## PC_ 3
## Positive: Fth1, Ubb, Ftl1, Fxyd2, Ldhb, Car15, Prdx1, Cd63, Eif1, Cd9
## Mt1, Mpc2, Mgst1, Tmem213, Clu, Bsg, Mdh1, Aldoa, Ramp3, Itm2b
## Selenow, S100a1, Tmem59, Spp1, Tspo, Tmem176b, Tle5, Atpif1, Wfdc2, Tmbim6
## Negative: Mir6236, Egf, CT010467.1, Umod, Tmem52b, Neat1, Nme7, Malat1, Slc12a1, Fos
## Kcnq1ot1, Lars2, Sult1d1, Jun, Foxq1, Etnk1, Dst, Wnk1, Slc5a3, Junb
## Atp1b1, Atrx, Pnisr, Egr1, Gls, Ivns1abp, Syne2, Hnrnpa2b1, Chka, Srrm2
## PC_ 4
## Positive: Pappa2, Aard, Tmem52b, Umod, Egf, Sult1d1, Tmem158, Foxq1, Mcub, Ptgs2
## Dctd, Wwc2, Iyd, Car15, Ramp3, Ptprz1, Hsp90b1, Cd9, S100g, Wnk1
## Pth1r, Cdkn1c, Defb42, Clu, Tmsb4x, Ctsc, Tagln2, Col4a3, Bmp2, Tsc22d1
## Negative: Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Ptger3, Neat1, Ckb
## Fgf9, Egfl6, Ivns1abp, Mfsd4a, Defb1, Car4, CT010467.1, Plet1, Fkbp11, Atp5md
## Chchd10, Atp5k, Igfbp5, Cox6c, Tmem213, Chka, Mgst3, Avpr1a, Atp5mpl, Atp1a1
## PC_ 5
## Positive: Hspa1a, Hspa1b, CT010467.1, Klk1, Pappa2, Car15, Cldn10, Fth1, Jun, Mir6236
## Lars2, Hspa8, Wfdc2, Itm2b, Fau, Ftl1, Uba52-ps, Pik3r1, Sfrp1, Ly6a
## Ptger3, Eef1a1, Mal, Atp1a1, Aard, Tspo, Tmem176a, Hsp90aa1, Cited2, Gpx4
## Negative: S100g, Actb, Tmem52b, Sdc4, Tmsb10, Egr1, Uroc1, Abhd2, Ppia, Serf2
## Atf3, Slc39a1, Alkbh5, Igfbp7, Ndufb1-ps, Atp5md, Cebpd, Cox6c, Ramp3, Ndufa3
## Gnb1, Atp5k, Atp5e, Wfdc15b, Kdm6b, Rbm47, Ldhb-ps, Atp5mpl, Fam107a, Ranbp3l## Computing nearest neighbor graph## Computing SNN## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11455
## Number of edges: 390053
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7071
## Number of communities: 27
## Elapsed time: 0 seconds## 21:39:15 UMAP embedding parameters a = 0.9922 b = 1.112## 21:39:15 Read 11455 rows and found 36 numeric columns## 21:39:15 Using Annoy for neighbor search, n_neighbors = 30## 21:39:15 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 21:39:15 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpumjlXP/file671a18ecb892
## 21:39:15 Searching Annoy index using 1 thread, search_k = 3000
## 21:39:16 Annoy recall = 100%
## 21:39:17 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 21:39:17 Initializing from normalized Laplacian + noise (using RSpectra)
## 21:39:18 Commencing optimization for 200 epochs, with 514278 positive edges
## 21:39:18 Using rng type: pcg
## 21:39:20 Optimization finished
DotPlot(SO2,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()## Warning: Could not find Slc5a12 in the default search locations, found in 'RNA'
## assay instead## Warning: Could not find Kdr in the default search locations, found in 'RNA'
## assay instead## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b
6 Finding Markers
Idents(SO4) <- "seurat_clusters"
all_markers <- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 0.5)## Calculating cluster 0## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session## Calculating cluster 1## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 5## Calculating cluster 6## Calculating cluster 7## Calculating cluster 8## Calculating cluster 9## Calculating cluster 10## Calculating cluster 11## Calculating cluster 12## Calculating cluster 13## Calculating cluster 14## Calculating cluster 15## Calculating cluster 16## Calculating cluster 17## Calculating cluster 18## Calculating cluster 19## Calculating cluster 20## Calculating cluster 21## Calculating cluster 22## Calculating cluster 23## Calculating cluster 24## Calculating cluster 25## Calculating cluster 26all_markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
slice_head(n = 3) %>%
ungroup() -> top3
DoHeatmap(SO4, features = top3$gene) + NoLegend()
7 Assigning Subclasses for clusters
SO4@meta.data <- SO4@meta.data %>%
mutate(
manual_groups = dplyr::case_when(
seurat_clusters %in% c(0, 1, 3, 5, 6, 12, 13, 16,19) ~ "type_1a",
seurat_clusters %in% c(2,24) ~ "type_1b",
seurat_clusters == 7 ~ "type_2a",
seurat_clusters == 9 ~ "type_3b",
seurat_clusters %in% c(10,18,4) ~ "type_1c",
seurat_clusters %in% c(11, 14,15,20,23,8) ~ "type_3a",
seurat_clusters %in% c(17,26) ~ "type_2b",
seurat_clusters == 21 ~ "type_4",
seurat_clusters == 22 ~ "type_5",
seurat_clusters == 25 ~ "type_6",
)
)
DimPlot(object = SO4, reduction = "umap", group.by = "manual_groups", label = TRUE)
Idents(SO4) <- "manual_groups"
manual_markers<- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 1)## Calculating cluster type_3a## Calculating cluster type_1a## Calculating cluster type_1b## Calculating cluster type_1c## Calculating cluster type_3b## Calculating cluster type_4## Calculating cluster type_2a## Calculating cluster type_2b## Calculating cluster type_6## Calculating cluster type_5manual_markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 0.5, p_val_adj < 0.05) %>%
slice_head(n = 3) %>%
ungroup() -> top3m
DoHeatmap(SO4, features = top3m$gene,group.by = "manual_groups") + NoLegend()
# Assign your markers as a vector
genes_to_plot <- c(
"Fos", "Socs3", "Gadd45b", "Wee1", "Hspa1a", "Sat1",
"Cxcl12", "Itpr2", "Bmp4", "Casr", "Grin2c", "Irx3", "Rap1gap", "App", "Wwc1",
"Syt5", "Syn3", "Cacna1d", "Slc6a7", "Robo2", "Begain",
"Pappa1", "Nos1", "Ptgs2", "Bmp2", "Atp2a3"
)
DoHeatmap(SO4, features = genes_to_plot, group.by = "manual_groups")## Warning in DoHeatmap(SO4, features = genes_to_plot, group.by =
## "manual_groups"): The following features were omitted as they were not found in
## the scale.data slot for the SCT assay: Pappa1, Wwc1, Irx3, Grin2c, Casr
Idents(SO4) <- "manual_groups"
DotPlot(SO4,
features = genes_to_plot,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()## Warning: The following requested variables were not found: Pappa1
# Supergroup
SO4@meta.data <- SO4@meta.data %>%
mutate(
supergroup = case_when(
manual_groups %in% c("type_1a", "type_1b", "type_1c") ~ "type_1",
manual_groups %in% c("type_2a", "type_2b") ~ "type_2",
manual_groups %in% c("type_3a", "type_3b") ~ "type_3",
manual_groups == "type_4" ~ "type_4",
manual_groups == "type_5" ~ "type_5",
manual_groups == "type_6" ~ "type_6"
)
)
DimPlot(object = SO4, reduction = "umap", group.by = "supergroup", label = TRUE)
Idents(SO4) <- "supergroup"
supergroup_markers<- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 1)## Calculating cluster type_3## Calculating cluster type_1## Calculating cluster type_4## Calculating cluster type_2## Calculating cluster type_6## Calculating cluster type_5supergroup_markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 0.5, p_val_adj < 0.05) %>%
slice_head(n = 5) %>%
ungroup() -> top5s
DoHeatmap(SO4, features = top5s$gene,group.by = "supergroup") + NoLegend()
# DotPlot supergroups
# Rename the supergroup (assuming the column is called 'supergroup')
SO4@meta.data$supergroup <-
factor(SO4@meta.data$supergroup,
levels = c("type_1", "type_2", "type_3", "type_4", "type_5","type_6"))
Idents(SO4) <- "supergroup"
markers.to.plot2 <- c("Slc12a1","Nos1","Pappa2","Fos","Egf","Cxcl10","Il1f6","Ifitm3")
DotPlot(SO4,
features = markers.to.plot2,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
8 DotPlot Manual Groups
markers.to.plot3 <- c("S100g","Ptger3","Pappa2","Hspa1b","Fos","Egf","Foxq1","Cxcl10","Il1f6","Ifitm3")
SO4@meta.data$manual_groups <-
factor(SO4@meta.data$manual_groups,
levels = c("type_1a","type_1b","type_1c","type_2a", "type_2b","type_3a", "type_3b", "type_4", "type_5","type_6"))
Idents(SO4) <- "manual_groups"
DotPlot(SO4,
features = markers.to.plot3,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
library(dplyr)
all_markers <- all_markers %>%
arrange(desc(avg_log2FC)) %>%
dplyr::select(gene, everything())
# Split by cluster (ident column)
marker_list <- split(all_markers, all_markers$cluster)
# Sort names alphanumerically
sorted_cluster_names <- names(marker_list)[order(as.numeric(names(marker_list)))]
# Create a workbook
wb <- createWorkbook()
# Add each cluster as a new worksheet
for (cluster_name in sorted_cluster_names) {
addWorksheet(wb, sheetName = cluster_name)
writeData(wb, sheet = cluster_name, marker_list[[cluster_name]])
}
date <- format(Sys.Date(), "%Y%m%d")
# Save workbook
saveWorkbook(wb, here("jk_code", paste0(date, "_", "_SO4_FindAllMarkers_By_Cluster2.xlsx")), overwrite = TRUE)
# Saving Seurat Object
save(SO4, file = here("jk_code", "JK_clean_MD_2.rds"))